def gmm_loglike(data,alpha,mu,sigma):
loglike_tmp=np.zeros((data.shape[0],alpha.shape[0]))
for j in range(alpha.shape[0]):
loglike_tmp[:,j]=alpha[j]*stats.multivariate_normal.pdf(data,mu[j],sigma[j])
loglike_tmp=np.apply_along_axis(lambda x: np.log(np.sum(x)),1,loglike_tmp)
loglike=np.sum(loglike_tmp)
return(loglike)
def Gaussian_Expect(data,alpha,mu,sigma, clusters=None,*,allowsingular:bool=False,normalize:bool=True,method='numerical'):
clusters=alpha.shape[0] if clusters is None else clusters
if method=='scipy':
tau=np.zeros((data.shape[0],clusters))
for j in range(clusters):
tau[:,j]=alpha[j]*stats.multivariate_normal.pdf(data,mu[j],sigma[j],allow_singular=bool(allowsingular))
if normalize:
normalizer=np.sum(tau,1)
tau=tau/normalizer.reshape(-1,1)
elif method=='numerical':
log_tautild=np.zeros((data.shape[0],clusters))
log_sumtautild=np.zeros((data.shape[0],clusters))
log_tau=np.zeros((data.shape[0],clusters))
for j in range(clusters):
#On ajoute pour chaque cluster J la log-vraisemblance de x.
#On va diagonaliser Sigma
vals, vecs = np.linalg.eigh(sigma[j])
sign, logdet = np.linalg.slogdet(sigma[j])
valsinv = 1./vals
U = vecs * np.sqrt(valsinv)
dev = data - mu[j]
#On passe par la distance de Mahalanobis
maha = np.sum(np.square(np.dot(dev, U)),1)
log_tautild[:,j]=np.log(alpha[j]) -0.5 * (maha + logdet)
log_sumtautild=np.log(np.sum(np.exp(log_tautild),1))
if normalize:
tau=np.exp(log_tautild)/np.exp(log_sumtautild).reshape(-1,1)
else:
tau=np.exp(log_tautild)
else:
raise ValueError('error, no methods chosen')
return(tau)
def Gaussian_Maxim(data,tau):
sigma=[]
mu=[]
alpha=np.sum(tau,axis=0)
alpha=alpha/np.sum(alpha)
for j in range(tau.shape[1]):
mu.append(np.apply_along_axis(lambda x: np.average(x,weights=tau[:,j]),0,data))
sigma.append((tau[:,j].reshape(-1,1)*(data-mu[j])).T@(data-mu[j])/np.sum(tau[:,j]))
return(dict(alpha=alpha,mu=mu,sigma=sigma))
def EM_Gaussian(data,max_iter=100,clusters=3,method='numerical'):
kmeans=cluster.KMeans(n_clusters=clusters,init='k-means++',n_init='auto').fit(data)
labels=kmeans.labels_
loglike_seq=np.zeros(max_iter+1)
mu=[]
sigma=[]
alpha=np.ones(clusters)
for i in range(clusters):
index=np.where(labels==i)[0]
mu.append(np.apply_along_axis(np.mean,0,data[index,:]))
sigma.append(np.cov([data[index,d] for d in range(data.shape[1])]))
alpha[i]=len(index)/data.shape[0]
for i in range(max_iter):
loglike_seq[i]=gmm_loglike(data,alpha,mu,sigma)
tau=Gaussian_Expect(data,alpha,mu,sigma,clusters,method=method)
results=Gaussian_Maxim(data,tau)
mu=results["mu"]
alpha=results["alpha"]
sigma=results["sigma"]
loglike_seq[-1]=gmm_loglike(data,alpha,mu,sigma)
clusterid=np.apply_along_axis(np.argmax,1,tau)
return(dict(alpha=alpha,mu=mu,sigma=sigma,tau=tau, loglike=loglike_seq,theta=[alpha,mu,sigma],clusterid=clusterid))
def Gaussian_contour_2d(data,results_EM,title="Clusters obtenus par algorithme EM",x=None,y=None,ax=None):
if ax is None:
fig, ax=plt.subplots(figsize=(15,15))
#plt.figure(figsize=(10,10))
ax.plot(data[:,0],data[:,1],'o',alpha=.6,lw=.8)
for i in range(results_EM['alpha'].shape[0]):
#On plotte uniquement les gaussiennes dans le cluster prédit pour éviter un comportement aberrant
index=np.where(results_EM['clusterid']==i)[0]
xmin = min(data[index,0])
xmax = max(data[index,0])
ymin = min(data[index,1])
ymax = max(data[index,1])
X=np.linspace(xmin,xmax)
Y=np.linspace(ymin,ymax)
X,Y=np.meshgrid(X,Y)
z=stats.multivariate_normal(results_EM['mu'][i][:2],results_EM['sigma'][i][:2,:2]).pdf(np.dstack((X,Y)))
ax.contour(X, Y, z,levels=9)
ax.set_title(title)
ax.set_xlabel(x)
ax.set_ylabel(y)
